#include "Function.hpp"
#include "EquationSolver.hpp"
#include <iostream>
#include <cmath>

const double Pi = acos(-1.); 

class F : public Function{
public: 
	double L = 10.;
	double r = 1.; 
	double V = 12.4;
	double operator() (double h) const{
		
		return L*(0.5*Pi*pow(r,2) - pow(r,2)*asin(h/r) - h*pow(r*r-h*h,0.5)) - V;
	}
	double derivative(double h) const{
		return L*(-r/(1 - pow(h/r,2)) - pow(r*r-h*h,0.5) - pow(h,2)*pow(r*r - h*h,-0.5));
	}
};  

void solvef_1(double a, double b, double eps){
	std::cout << "solve by bisection" << "from ["<< a <<","<< b<<"]"<< std::endl;
	Bisection_Method solver1(F(), a, b);
	double x = solver1.solve(); 
	std::cout << "the root is " << x <<  std::endl;	
}

void solvef_2(double x0, double eps){
  std::cout << "solve by newton near "<< x0 << std::endl;
	Newton_Method solver2(F(), x0, eps, 1000);
	double x = solver2.solve(); 
	std::cout << "the root is " << x << std::endl;
}

void solvef_3(double x0, double x1, double eps){
  std::cout << "solve by secant with initial"<< x0 <<"and"<< x1 << std::endl;
	Secant_Method solver3(F(), x0, x1, eps);
	double x = solver3.solve(); 
	std::cout << "the root is " << x << std::endl;
}

int main(){
	double eps = 0.01;
	solvef_1(0, 1, eps);
	solvef_2(0.9, eps);
	solvef_3(0.9, 0.1, eps); 
}
